clc; clear; close all;

% 参数设置
g11 = 1.4; g12 = -2; g21 = 1.4; g22 = 1.6;

k_list = linspace(0, 3, 900);

dt = 0.01;
t_total = 400;
t_transient = 350;

% 两组初始条件
IC1 = [0.01; 0; 0];
IC2 = [-0.01; 0; 0];

% 用cell存储，方便parfor
peaks_IC1 = cell(length(k_list), 1);
kvals_IC1 = cell(length(k_list), 1);

peaks_IC2 = cell(length(k_list), 1);
kvals_IC2 = cell(length(k_list), 1);

parfor idx = 1:length(k_list)
    k = k_list(idx);

    % IC1积分
    [T1, X1] = ode45(@(t,x)MBNHNN(x,k,g11,g12,g21,g22), [0 t_total], IC1);
    start_idx1 = find(T1 >= t_transient, 1, 'first');
    x1_steady1 = X1(start_idx1:end, 1);
    [pks1, ~] = findpeaks(x1_steady1);
    peaks_IC1{idx} = pks1;
    kvals_IC1{idx} = repmat(k, length(pks1), 1);

    % IC2积分
    [T2, X2] = ode45(@(t,x)MBNHNN(x,k,g11,g12,g21,g22), [0 t_total], IC2);
    start_idx2 = find(T2 >= t_transient, 1, 'first');
    x1_steady2 = X2(start_idx2:end, 1);
    [pks2, ~] = findpeaks(x1_steady2);
    peaks_IC2{idx} = pks2;
    kvals_IC2{idx} = repmat(k, length(pks2), 1);
end

% 合并结果
x1_peaks_IC1 = vertcat(peaks_IC1{:});
k_vals_IC1 = vertcat(kvals_IC1{:});
x1_peaks_IC2 = vertcat(peaks_IC2{:});
k_vals_IC2 = vertcat(kvals_IC2{:});

% 绘图
figure;
hold on;
plot(k_vals_IC1, x1_peaks_IC1, '.', 'Color', [0.85 0.33 0.1], 'DisplayName', 'IC_1 = (0.01,0,0)');
plot(k_vals_IC2, x1_peaks_IC2, '.', 'Color', [0 0.45 0.74], 'DisplayName', 'IC_2 = (-0.01,0,0)');
xlabel('k', 'FontSize', 14);
ylabel('x_{1peak}', 'FontSize', 14);
grid on;
hold off;


